Technical  Report  HL-95-9 
September  1995 


WwW 


US  Army  Corps 
of  Engineers 

Waterways  Experiment 
Station 


Preliminary  Deveiopment  of  a  Three- 
Dimensionai  Numericai  Modei  for 
Reservoir  Hydrodynamics 

by  Robert  S.  Bernard 

I 


Approved  For  Public  Release;  Distribution  Is  Unlimited 


19960124  116 

DTIC  QU/Om  IKSPECJTED  1 


Prepared  for  U.S.  Army  Engineer  District,  Chicago 


The  contents  of  this  report  are  not  to  be  used  for  advertising, 
publication,  or  promotional  purposes.  Citation  of  trade  names 
does  not  constitute  an  official  endorsement  or  approval  of  the  use 
of  such  conmiercial  products. 


PRINTED  ON  RECYCLED  PAPER 


Technical  Report  HL-95-9 
September  1 995 


Preliminary  Development  of  a  Three- 
Dimensional  Numerical  Model  for 
Reservoir  Hydrodynamics 

by  Robert  S.  Bernard 

U.S.  Army  Corps  of  Engineers 
Waterways  Experiment  Station 
3909  Halls  Ferry  Road 
Vicksburg,  MS  39180-6199 


Final  report 

Approved  for  public  release;  distribution  is  unlimited 


Prepared  for  U.S.  Army  Engineer  District,  Chicago 
Chicago,  IL  60606-7206 


Waterways  Experiment  Station  Cataloging-in-Publication  Data 

Bernard,  Robert  S. 

Preliminary  development  of  a  three-dimensional  numerical  model  for  reservoir 
hydrodynamics  /  by  Robert  S.  Bernard  ;  prepared  for  U.S.  Army  Engineer  District, 
Chicago. 

92  p. ;  ill. ;  28  cm.  —  (Technical  report ;  HL-95-9) 

1 .  Reservoirs — Mathematical  models.  2.  Hydrodynamics — Mathematical  models. 
3.  Turbulence  —  Mathematical  models.  I.  U.S.  Army.  Corps  of  Engineers.  Chicago 
District.  II.  U.S.  Army  Engineer  Waterways  Experiment  Station.  III.  Hydraulics 
Laboratory  (U.S.  Army  Engineer  Waterways  Experiment  Station)  IV.  Title.  V.  Series: 
Technical  report  (U.S.  Army  Engineer  Watenways  Experiment  Station) ;  HL-95-9. 
TA7  W34  no.HL-95-9 


Contents 


Preface  .  v 

Conversion  Factors,  Non-SI  to  SI  Units  of  Measurement . vi 

1—  Introduction  .  1 

Background .  1 

Purpose  and  Scope  .  3 

2—  Governing  Equations .  4 

Incompressible  Flow .  4 

Passive  Transport .  7 

Heat  Transfer  .  8 

Turbulence .  9 

Curvilinear  Coordinates  . 12 

3—  Numerical  Considerations . 15 

Finite-Volume  Discretization . 15 

Computational  Coordinates  . 16 

Face-Centered  Fluxes  . 17 

Cartesian  Velocities . 17 

Predictor-Corrector  Scheme  for  Updating  Fluxes . 18 

Predictor-Corrector  Scheme  for  Transport  Equations . 19 

Single-Step  Upwind  Scheme  for  Turbulence . 21 

Inflow/Outflow  Boundaries . 22 

4—  Model  Operation  . 24 

Composite  Grids . 24 

Grid  Generation . 24 

Inflow/Outflow  Data . 25 

Initial  Conditions  . 25 

Flow  Development  . 25 

Flow  Visualization  . 26 


iii 


5 —  Flow  Simulations . 27 

Provisional  Grid  for  McCook  Reservoir . .  .  27 

Imposed  Flow  Conditions . 28 

Presentation  of  Computed  Results . 29 

Cold  Reservoir  with  Cold  Inflow  . 29 

Warm  Reservoir  with  Cold  Inflow . 29 

Stratified  Reservoir  with  Cold  Inflow . 30 

6—  Conclusions  and  Recommendations . . . 31 

References  . 33 

Figures  1-50 

Appendix  A;  Notation .  A.1 

SF  298 


IV 


Preface 


This  report  presents  the  results  of  research  and  development  conducted 
from  April  1991  through  May  1994  by  personnel  of  the  Reservoir  Water 
Quality  Branch  (RWQB),  Hydraulic  Structures  Division  (HSD),  Hydraulics 
Laboratory  (HL),  U.S.  Army  Engineer  Waterways  Experiment  Station  (WES). 
Funding  was  provided  by  the  U.S.  Army  Engineer  District,  Chicago. 

Dr.  Robert  S.  Bernard,  RWQB,  prepared  this  report  under  the  general 
supervision  of  Messrs.  Frank  A.  Herrmann,  Jr.,  Director,  HL;  Richard  A. 
Sager,  Assistant  Director,  HL;  and  Glenn  A.  Pickering,  Chief,  HSD. 
Technical  counsel  was  provided  by  Messrs.  Mike  Schneider,  Stacy 
Howington,  and  Steve  Wilhelms,  RWQB;  by  Drs.  Hartmut  Kapitza  and  Dieter 
Eppel,  GKSS  Research  Centre,  Geesthacht,  Germany;  and  by  Drs.  Abi 
Arabshahi,  Lafe  Taylor,  Mike  Stokes,  Steve  Davis,  Boyd  Gatlin,  Dave 
Whitfield,  Joe  Thompson,  Roger  Briley,  Pasquale  Cinnella,  Dave  Huddleston, 
Bharat  Soni,  and  Tim  Swafford,  and  Mr.  Mike  Remotigue,  Engineering 
Research  Center  for  Computational  Field  Simulation,  Mississippi  State 
University.  Dr.  Kapitza  provided  the  conjugate-gradient  solver  used  for  the 
pressure  Poisson  equation. 

Technical  monitors  for  the  Chicago  District  were  Mr.  Tom  Fogarty  and 
Mrs.  Linda  Sorn. 

At  the  time  of  preparation  of  this  report.  Dr.  Robert  W.  Whalin  was 
Director  of  WES,  and  COL  Bruce  K.  Howard,  EN,  was  Commander. 


The  contents  of  this  report  are  not  to  be  used  for  advertising,  publication, 
or  promotional  purposes.  Citation  of  trade  names  does  not  constitute  an 
official  endorsement  or  approval  of  the  use  of  such  commercial  products. 


Conversion  Factors, 
Non-SI  to  SI  Units  of 
Measurement 


Non-SI  units  of  measurement  used  in  this  report  can  be  converted  to  SI  units 
as  follows; 


Multiply 

By 

To  Obtain 

cubic  feet 

0.02831685 

cubic  meters 

degrees  (angle) 

0.01745329 

radians 

feet 

0.3048 

meters 

square  feet 

0.09290304 

square  meters 

1  Introduction 


Background 


Computational  fluid  dynamics  (CFD)  is  now  used  routinely  by  researchers 
and  design  engineers  alike.  New  CFD  codes  (numerical  flow  models)  appear 
every  year.  Although  some  of  these  are  applicable  for  broad  classes  of  flow 
problems,  no  single  model  is  likely  to  gain  universal  acceptance  in  the 
foreseeable  future.  No  one  has  yet  written  a  general-purpose  CFD  code  that 
works  reliably  and  efficiently  at  all  Reynolds  numbers,  Froude  numbers,  and 
Mach  numbers.  In  spite  of  monumental  gains  in  computer  power,  numerical 
flow  modeling  remains  a  problem-specific  art. 

Buoyant  flow  in  reservoirs  falls  in  a  class  of  problems  for  which  the  Mach 
number  is  very  small.  Thus,  the  density  can  be  assumed  to  vary  only  slightly 
with  temperature  and  not  at  all  with  pressure.  The  resulting  density  gradients 
are  so  small  that  the  governing  equations  reduce  to  the  incompressible  Navier- 
Stokes  equations  with  a  vertical  perturbation  added  for  buoyancy.  The  Froude 
number  is  generally  small  enough  in  reservoirs  to  justify  neglect  of  surface 
waves,  and  the  free  surface  can  be  assigned  a  uniform  vertical  velocity.  In 
contrast,  the  Reynolds  number  may  exceed  10^,  and  empirical  corrections  are 
needed  to  model  the  small-scale  influence  of  turbulence. 

Some  buoyant  flow  models  assume  the  pressure  in  a  water  colunm  to  be 
hydrostatic,  i.e.,  proportional  to  the  weight  of  the  water  itself.  This 
eliminates  the  vertical  component  of  the  momentum  equation  and  guarantees 
that  the  remaining  equations  will  all  be  of  the  same  mathematical  type,  but  it 
is  acceptable  only  when  the  vertical  acceleration  is  negligible.  Deep 
reservoirs  with  submerged  structures  and  sharply  varying  topography 
generally  render  the  hydrostatic  assumption  invalid. 

MAC3D  is  a  numerical  model  for  buoyant  incompressible  flow,  with 
emphasis  on  reservoir  hydrodynamics.  In  general,  the  model  is  intended  for 
hydraulic  applications  in  deep  water  where  vertical  motion  and  vertical 
acceleration  are  both  important.  It  uses  a  variant  of  MacCormack’s  method 
(MacCormack  1969;  Bernard  1992)  to  solve  the  Reynolds-averaged  Navier- 
Stokes  equations  for  three-dimensional  (3-D)  incompressible  flow,  which  are 
discretized  with  six-sided  finite-volume  cells.  The  grid  cells  can  be 
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nonrectangul'ar  in  cartesian  {x,y,z)  space,  but  they  are  uniformly  rectangular  in 
the  computational  (ij.k)  space  where  they  are  indexed  by  the  integer 
coordinates  (iJ.k).  This  sequential  ordering  of  grid  cells  in  an  ijk-aii&y  is 
commonly  called  a  structured  grid. 

To  improve  local  resolution  and  reduce  computer  memory  requirements, 
MAC3D  uses  composite  (multiblock)  structured  grids,  i.e.,  structured  grids 
composed  of  one  or  more  subdomains  (also  known  as  blocks).  The  individual 
subdomains  are  rectangular  in  the  computational  (i,j,k)  space,  but  they  can 
have  any  shape  and  orientation  in  cartesian  ix,y,z)  space  as  long  as  their  grid 
nodes  match  along  the  shared  boundaries  between  the  subdomains.  The 
discrete  flow  variables  are  staggered,  with  contravariant  (normal)  velocity 
components  defined  at  the  centroids  of  the  cell  faces,  and  scalar  quantities 
(pressure,  temperature,  etc.)  defined  at  the  cell  centers. 

MAC3D  accounts  for  turbulence  by  using  a  k-e  turbulence  model  (Launder 
and  Spalding  1974),  which  consists  of  two  semiempirical  equations  for  the 
production  and  transport  of  the  turbulence  energy  k  and  the  turbulence 
dissipation  rate  €.  Buoyancy  arises  from  temperature-dependent  density 
variations  that  add  a  perturbing  force  to  the  vertical  momentum  equation  in 
the  presence  of  gravity.  The  transport  of  heat  (temperature)  and  other  passive 
constituents  is  governed  by  advection-diffusion  equations,  in  which  the  dif¬ 
fusion  coefficients  are  proportional  to  the  eddy  viscosity  obtained  from  the 
turbulence  model. 

/ 

The  normal  velocity  along  inflow/outflow  boundaries  can  be  imposed  either 
as  an  unchanging,  user-specified  distribution  or  as  a  time- varying  distribution 
extrapolated  from  the  flow  just  inside  the  grid.  In  the  latter  case,  the 
extrapolated  velocities  are  obtained  from  a  discrete  radiation  condition 
proposed  by  Orlanski  (1976),  which  transmits  internal  disturbances  out  of  the 
grid  with  negligible  reflection.  Solid  boundaries  can  be  designated  as  either 
no-slip  (frictional)  or  slip  (frictionless).  The  shear  stress  along  no-slip 
boundaries  is  obtained  either  from  Manning’s  equation  or  from  an  empirical 
wall  function  that  defines  the  relation  between  velocity  and  distance  from  the 
boundary.  Free  surfaces  are  idealized  as  (rigid)  slip  boundaries. 

MAC3D  tentatively  assumes  all  solid  boundaries  and  free  surfaces  to  be 
impermeable  and  adiabatic  (thermally  insulated),  so  that  heat  and  other  passive 
constituents  can  enter  or  leave  the  flow  field  only  through  the  inflow/outflow 
boundaries.  The  actual  rate  of  heat  transfer  through  a  free  surface  depends  on 
external  factors  such  as  radiation  and  wind  (which  also  creates  momentum 
transfer),  and  these  will  be  dealt  with  in  future  work.  For  the  time  being, 
however,  the  model  allows  the  adiabatic  constraint  to  be  relaxed  in  the 
following  manner. 

Suppose  that  a  reservoir  is  initially  stratified  with  cold  water  at  the  bottom 
and  warm  water  at  the  surface.  If  all  the  boundaries  are  adiabatic  (no  heat 
transfer)  and  impermeable  (no  inflow  or  outflow),  then  diffusion  (heat  conduc¬ 
tion)  will  gradually  eliminate  the  initial  stratification,  leaving  the  reservoir 
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with  a  single  unifonn  temperature.  In  some  cases,  however,  one  may  wish  to 
know  only  the  manner  in  which  a  particular  flow  alters  the  initial  stratifica¬ 
tion,  and  not  the  rate  at  which  the  insulated  reservoir  otherwise  comes  to 
equilibrium.  Accordingly,  if  the  user  so  chooses,  MAC3D  will  ignore  any 
change  in  temperature  not  produced  by  the  flow  itself.  In  other  words,  it  will 
ignore  background  heat  conduction.  This  allows  the  user  to  isolate  the  effect 
of  the  flow  upon  the  stratification.  Otherwise,  the  model  will  allow  the 
temperature  to  evolve  toward  an  equilibrium  state  dictated  by  the  imposed 
flow,  the  full  heat  conduction,  and  the  adiabatic  boundaries. 


Purpose  and  Scope 


The  planned  McCook  Reservoir  near  Chicago,  IL,  will  serve  as  a  holding 
tank  for  runoff  and  sewage  (prior  to  treatment)  through  the  middle  of  the 
twenty-first  century.  Its  design  and  construction  pose  major  engineering  chal¬ 
lenges  because  of  its  size— roughly  4,700  ft^  long,  3,300  ft  wide,  and  110  to 
150  ft  deep.  Safe  operation  of  the  reservoir  demands  an  intimate  understand¬ 
ing  of  its  hydrodynamics,  and  this  warrants  the  development  of  a  nonhydro¬ 
static,  3-D  numerical  model  as  a  means  to  that  end.  A  model  of  this  kind 
would  be  useful  not  only  for  McCook  Reservoir,  but  for  other  applications  as 
well. 

This  report  documents  the  preliminary  development  of  the  MAC3D 
numerical  model  and  offers  qualitative  demonstrations  of  MAC3D’s 
capabilities  using  a  provisional  10-block  grid  for  McCook  Reservoir. 
Quantitative  testing  and  empirical  validation  of  the  model  are  already  under¬ 
way,  but  these  are  lengthy  efforts  that  will  be  documented  in  future  reports. 

Chapter  2  outlines  the  governing  equations  for  the  major  flow  variables, 
while  Chapter  3  describes  the  discrete  methods  used  for  their  numerical  solu¬ 
tion.  Chapter  4  summarizes  briefly  the  operation  of  the  model  as  a  whole, 
and  Chapter  5  presents  the  results  of  preliminary  flow  calculations  for 
McCook  Reservoir.  Chapter  6  offers  conclusions  and  recommendations  with 
regard  to  further  model  development. 


*  A  table  of  factors  for  converting  non-SI  units  of  measurement  to  SI  units  is  found  on 
page  vi. 
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2  Governing  Equations 


Incompressible  Flow 

The  speed  of  sound  in  fresh  water  is  about  4,700  ft/sec.  With  Mach 
number  defined  as  the  ratio  of  flow  speed  to  sound  speed,  water  flowing  at 
47  ft/sec  has  a  Mach  number  of  0.01 .  At  Mach  numbers  on  the  order  of 
10'^,  water  behaves  as  though  it  were  truly  incompressible  (except  for 
questions  related  to  acoustics).  In  this  context,  an  incompressible  fluid  is 
defined  as  one  whose  density  varies  with  temperature  but  not  with  pressure. 

The  density  of  liquid  water  is  greatest  at  4°  C,  above  which  it  falls 
monotonically  with  temperature  up  to  100°  C.  This  behavior  of  the  density  p 
is  conveniently  expressed  as 


p  =  po[i  -  e{T)] 


(1) 


where  ^ 

Pq  =  reference  (maximum)  density 
6  =  relative  deviation  from  the  reference  value 
T  =  temperature 

For  temperatures  in  the  range  0-30°  C,  B  can  be  approximated  with 
e  =A\T-To\  +  B{T-Tq)^ 
where 

A  =  5.6250  X  10-^ 

^  For  convenience,  symbols  and  abbreviations  are  listed  in  the  notation  (Appendix  A). 
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At  30  °C,  Equation  2  gives  6  =  0.00458,  indicating  that  the  density 
deviates  only  slightly  from  its  reference  value.  Dissolved  constituents  such  as 
salt  can  also  produce  comparable  density  variations,  and  salinity  may  be 
included  in  future  versions  of  MAC3D.  In  any  case,  except  for  the  vertical 
buoyant  force  in  the  presence  of  gravity,  one  can  ignore  small  variations  of 
the  density  and  impose  the  fully  incompressible  constraint  for  conservation  of 
mass. 


div  u  -  0 


(3) 


where  the  symbol  div  denotes  the  divergence  operator,  u  represents  the 
velocity,  and  a  single  underbar  henceforth  indicates  a  vector.  Applying  the 
same  arguments  to  the  equation  for  conservation  of  momentum,  one  obtains 


Po^  =  divT  -  gradp  +  pg 


(4) 


where 
t  =  time 
T  =  shear  stress 
p  =  pressure 

g  =  acceleration  due  to  gravity 

and  a  double  underbar  henceforth  indicates  a  tensor.  The  symbol  grad 
represents  the  gradient  operator,  and  D/Dt  denotes  the  substantive  (total) 
time-derivative  operator, 

—  =  —  +  u-grad  (5) 

Dt  dt  - 


The  grad  and  div  operators  are  standard  vector  notation.  In  cartesian  (x.y.z) 
coordinates  they  take  the  following  form: 
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(6) 


..  du 
divu  =  —  + 
dx 


dv 


dw 

dz 


grad^  = 


d\f>  dip  3^ 
dx'  dy'  dz 


(7) 


The  triad  (m,v,>v)  gives  the  cartesian  components  of  the  vector  u,  and  the 
symbol  ^  here  denotes  any  scalar  fiinction  or  vector  component.  For  the 
latter,  the  operator  u- grad  yields 


, ,  d^  d\p 
u-gradxp  =  u-^  +  v-^ 
dx  dy 


+  w 


3^ 

dz 


(8) 


For  turbulent  flow,  the  velocity  in  Equations  3  and  4  is  the  Reynolds- 
averaged  velocity,  and  the  shear  stress  includes  the  viscous  stress  and  the 
turbulent  Reynolds  stress.  Invoking  the  Boussinesq  hypothesis  for  turbulence, 
the  two  stresses  are  combined  in  a  single  expression. 


=  M 


3m,-  dUj 

_ 1  +  _ i 

dxj  3x,- 


(9) 


The  dynamic  viscosity  fi  is  the  sum  of  the  molecular  viscosity  and  the 
eddy  viscosity  ixp  Since  assumes  in  practice  that  n  «  (ij. 

In  Equation  9,  the  subscripts  i  and  j  take  integer  values  1,  2,  and  3.  Thus, 
the  vector  (m^,M2,Mj)  represents  velocity  {u,v,w),  and  the  vector  (Xj,X2,x^) 
represents  position  (x,y,z)  in  a  cartesian  coordinate  system.  The  tensor  stress 
component  is  the  shear  stress  in  the  x/*  direction  on  a  plane  normal  to  the 
Xj^  direction  (and  vice  versa,  because  of  symmetry).  For  example,  tj2  is  the 
shear  stress  in  the  x-direction  on  a  plane  normal  to  the  y-direction,  given  by 


^12 


^21 


3mi 

^  3M2 

3m 

3^2 

3xj 

dy 

dx 

(10) 


Assuming  that  the  temperature,  viscosity,  and  pressure  are  known  at  each 
instant,  then  Equations  3  and  4  uniquely  determine  the  evolution  of  the  flow 
with  time.  Temperature  can  be  obtained  from  a  transport  equation  for  heat, 
and  eddy  viscosity  from  a  suitable  turbulence  model,  but  there  is  no  equation 
of  state  from  which  to  derive  the  pressure.  The  correct  pressure  gradient  for 
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incompressible  flow  is  that  which  ensures  conservation  of  mass  as  expressed 
by  Equation  3.  To  obtain  a  governing  equation  of  convenient  form  for  pres¬ 
sure,  one  takes  the  divergence  of  Equation  4.  This  produces  a  Poisson 
equation,  i.e., 


-div 

dx^  dy^  dz^ 


divT  -  pg 


(11) 


In  practice,  the  right-hand  side  of  Equation  1 1  is  obtained  in  the  following 
manner.  Suppose  that  the  flow  is  to  be  advanced  from  time  r  to  t  +  At.  Let 
u  be  the  mass-conserving  velocity  at  time  t,  and  let  m'  be  a  non-mass- 
conserving  velocity  computed  for  time  r  -I-  At  by  omitting  grad  p  from 
Equation  4.  Thus,  div  u'  ^  0,  where 


u 


—  {pQU'gradu  -  divr  -  pg] 
Po  ~  ~ 


(12) 


If  grad  p  is  the  (unknown)  pressure  gradient  that  has  to  be  imposed  at  time  t 
to  guarantee  conservation  of  mass  at  time  t  -I-  At,  then  Equation  1 1  reduces  to 


dz^  ^ 


In  other  words,  even  if  the  velocity  u  is  known  at  a  given  time,  the  pressure 
for  that  time  is  not  necessarily  known.  To  find  the  correct  pressure,  the 
momentum  equation  (4)  must  first  be  advanced  (integrated)  provisionally  over 
a  small  time  increment  At  without  including  a  pressure  gradient.  The  resulting 
mass-conservation  error  div  u'  ^  0  then  provides  the  right  side  of  the  Poisson 
equation  (11,  13)  whose  solution  yields  the  unknown  pressure  gradient. 


Passive  Transport 

The  transport  equations  for  quantities  other  than  momentum  all  have  the 
general  form  of  an  advection-diffusion  equation  with  a  source/sink  term  on  the 
right  side,  i.e.. 


D^P 

Dt 


=  divi^v^  grad  ^1/)  + 


(14) 
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where  is  the  diffiisivity  and  is  the  source/sink  term  for  the  transported 
quantity  i/.  One  ordinarily  assumes  to  be  proportional  to  the  kinematic 
fluid  viscosity  v,  i.e., 

V,  (15) 


where 

V  =  }!l  ^  ^  (16) 

P  PQ 

where  is  the  ratio  of  diffusivity  to  kinematic  viscosity  for  the  particular 
quantity  ip. 


Heat  Transfer 

For  incompressible  flow,  thermal  energy  (heat)  is  uncoupled  from  kinetic 
energy  and  is  directly  proportional  to  temperature.  Thus,  the  equation  for 
passive  transport  of  heat  reduces  to  an  equation  for  temperature,  i.e.. 


DT 

TJ 


=  div{Pjv  grad  T) 


+  S7 


(17) 


where  the  coefficient  Pj  is  the  ratio  of  thermal  diffusivity  to  kinematic 
viscosity,  also  known  as  the  Prandtl  number.  The  thermal  source/sink  term 
Sf  is  zero  in  the  absence  of  chemical,  biological,  or  radiative  processes  that 
generate  or  absorb  heat. 

Equation  17  determines  the  rate  of  internal  heat  transfer  through  the  fluid 
itself,  but  not  through  the  boundaries.  In  practice,  solid  boundaries  can  be 
regarded  as  adiabatic  walls  that  insulate  the  flow  from  its  surroundings.  The 
adiabatic  boundary  condition  for  temperature  is 


n-gradT  =  0 


(18) 


where  n  is  a  unit  vector  normal  to  the  boundary.  Equation  18  is  an  acceptable 
approximation  on  the  bottom  and  sides  of  a  reservoir,  but  generally  not  on  a 
free  surface.  Heat  transfer  across  a  free  surface  entails  radiation,  convection, 
and  conduction,  which  depend  on  sunlight,  wind,  and  ambient  air  tempera¬ 
ture.  When  the  influence  of  these  external  factors  is  unknown  or  poorly 
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defined,  the  following  procedure  offers  a  means  of  avoiding  some  of  the 
difficulties  arising  from  uncertainty. 

Suppose  a  reservoir  is  stably  stratified  initially  with  a  vertical  (z-dependent) 
temperature  distribution  T^(z)  that  is  constant  (or  very  slowly  changing)  with 
time.  One  may  think  of  Ts(z)  as  a  static  (background)  distribution  sustained 
by  unknown  conditions  outside  the  reservoir.  In  designing  artificial  measures 
to  overcome  or  eliminate  the  stratification,  it  is  expedient  to  compute 
immediate  response  to  a  perturbing  flow  and  to  ignore  slow  response  to 
external  conditions.  To  this  end.  Equations  17  and  18  are  replaced  by 


DT 

Wt 


div^PjV  grad  {T  -  7^)1 


(19) 


and 


n-grad(T  -  =  0 


(20) 


Equations  19  and  20  ensure  that  there  will  be  no  departure  from  the  static 
temperature  T^(z}  if  there  is  no  flow  (m  s  0)  in  the  reservoir.  This 
reformulation  of  the  equations  ignores  background  heat  conduction  (diffusion) 
and  seeks  to  isolate  the  influence  of  the  flow  in  particular  upon  stratification. 

External  influences  obviously  do  have  to  be  included  in  long-term  reservoir 
simulations,  but  these  have  been  deferred  for  later  work.  In  the  preliminary 
version  of  MAC3D,  Equation  19  is  used  by  default,  with  Equation  20 
imposed  on  all  boundaries  except  inflow/outflow  boundaries.  Background 
heat  conduction  (Equations  17  and  18)  is  optional. 


Turbulence 


Turbulence  arises  whenever  there  is  too  little  molecular  viscosity  to  prevent 
small  disturbances  from  growing  and  disrupting  a  laminar  flow.  The  process 
is  self-limiting,  so  that  a  turbulent  flow  can  be  thought  of  as  having  a  slowly 
varying  mean-flow  component  and  a  smaller,  rapidly  varying  turbulence 
component.  The  coupling  between  mean  flow  and  turbulence  is  so  strong  that 
one  cannot  be  calculated  independently  of  the  other,  but  the  resulting  range  of 
eddy  sizes  is  so  great  that  conventional  numerical  methods  cannot  resolve 
them  all.  Thus,  one  must  introduce  empirical  or  semiempirical  equations  to 
model  the  influence  of  turbulence  on  mean  flow,  and  vice  versa.  Henceforth, 
except  for  k  and  e,  all  reference  to  the  flow  variables  will  tacitly  imply  the 
mean  flow.  The  influence  of  buoyancy  in  the  production  and  dissipation  of 
turbulence  is  omitted  in  the  preliminary  version  of  MAC3D,  but  it  will  be 
incorporated  in  future  work. 
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One  of  the  most  widely  used  turbulence  models  is  the  A:-e  model  of 
Launder  and  Spalding  (1974).  Here  k  is  the  turbulence  energy  (per  unit 
mass),  and  e  is  its  dissipation  rate.  The  kinematic  eddy  viscosity  v  is  related 
to  k  and  e  by 


where  is  an  empirical  coefficient.  In  the  standard  k-e  model,  the  governing 
equations  for  k  and  e  are,  respectively, 

BJL  =  yF  -  e  +  —div(v gradk)  (22) 

Dt 

2 

—  =  C.vVi  -  Co—  +  —div(v grade)  (23) 

Dt  ^  k  ^  k 


The  first  term  I'F  on  the  right  side  of  Equation  22  is  the  production  rate  for 
turbulence  energy,  in  which 


T.  -1  2  2  2 

r  =  2\u^+Vy 


z  ^  {^z 


and  the  subscripts  x,  y,  and  z  indicate  partial  derivatives.  The  standard  set  of 
nondimensional  empirical  coefficients  is 

=  0.09 
Cj  =  1.44 
C2  =  1.92 
cr,  =  1.0 
ffj  =  1.3 

In  regions  of  low  velocity  and  high  turbulence  energy,  the  standard  k-e 
model  overpredicts  the  eddy  viscosity,  which  means  that  k  is  too  large,  or  e 
is  too  small,  or  both.  There  is  no  perfect  remedy  for  this  tendency,  but  ad  hoc 
modifications  that  preferentially  reduce  k,  enlarge  e,  or  damp  production  in 
low-velocity,  high-turbulence  regions  can  sometimes  be  tuned  to  improve  flow 
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predictions  in  general  (Bernard  1991,  1993).  In  MAC3D  the  quantity  T  is 
replaced  in  Equations  22  and  23  by  F*,  where 


r  *  =  r  tanh 


R*  v^|f 


(25) 


and  »  5  is  a  new  empirical  coefficient  whose  value  has  been  chosen  (by 
trial  and  error)  to  improve  predictions  for  two-dimensional  (2-D)  recirculating 
flow. 

The  turbulence  model  provides  equations  for  the  production  and  transport 
of  k  and  e,  but  not  for  the  necessary  boundary  conditions.  These  have  to  be 
deduced  from  other  considerations.  In  one  scenario  used  in  previous  work 
(Bernard  1993)  and  also  in  MAC3D,  the  condition  n*grarf  k  =  n  'grad  e  = 

0  is  imposed  on  all  impermeable  boundaries  except  no-slip  walls,  where 
tangential  velocity  and  turbulence  energy  are  both  assumed  proportional  to  the 
1/7  power  of  distance  from  the  wall.  The  tangential  shear  stress  on  no-slip 
walls  is  then  taken  to  be 


_  2^  «w 


(26) 


where  is  the  average  tangential  velocity  over  a  distance  5  normal  to  the 
wall.  The  quantity  is  a  wall  coefficient  whose  default  value  (unity)  applies 
only  for  hydraulically  smooth  walls,  but  this  value  can  be  changed  at  the 
user’s  discretion. 

In  an  alternative  scenario  offered  only  in  MAC3D,  the  condition  dk/dn  = 
de/dn  =  0  is  imposed  on  all  impermeable  boundaries,  and  the  shear  stress  on 
no-slip  walls  is  taken  to  be 

=  Po 


where  is  defined  as  before,  and  is  a  nondimensional  friction  coefficient. 
The  latter  is  obtained  from  the  following  relation  derived  from  Manning’s 
equation: 
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(28) 


14.5 


where  n  is  the  Manning  coefficient,  and  5  is  defined  as  before.  In  this  context 
n  serves  as  a  convenient  (and  familiar)  parameter  that  characterizes  the 
roughness  of  the  wall.  The  factor  14.5  applies  only  when  6  is  given  in  in 
Equation  28.  Otherwise,  it  becomes  9.8  when  b  is  given  in  meters. 

In  the  first  scenario,  the  energy  production  rate  pF*  is  calculated 
everywhere  using  the  definitions  for  F  and  F*  given  by  Equations  24  and  25. 
In  the  second  scenario,  however,  the  product  j'F*  is  replaced  adjacent  to  no¬ 
slip  boundaries  by 


3 

_  2  ^3/2 

* - ~T 

K  0 


(29) 


where  k  «  0.418  is  von  Karman’s  constant.  This  gives  the  production  rate  a 
value  that  it  would  otherwise  have  only  at  equilibrium  (when  production 
balances  dissipation  and  diffusion). 


where  the  subscripts  jc,y,z  and  indicate  partial  derivatives.  Using  the 
Gauss  theorem,  one  finds  that 
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div  M  =  A  (  ) 


(33) 


and  from  Equations  30-32,  it  follows  that 


u-gradyp  =  +  Wi^^)  (34) 


where  J  is  the  Jacobian  determinant, 


J  -  -  yj  y,  yf 

S(£,n,D  '  ’ 

h  h 


and  U,  V,  and  W  are  the  volumetric  fluxes  in  the  ^-,  rj-,  and  f-directions, 
respectively,  given  by 


U  =  +  ^yV  + 

(36) 

V  =  J(^ri^u  +  TiyV  + 

(37) 

(38) 

Lastly,  the  derivatives  of  the  curvilinear  coordinates  (^,?7,f)  with  respect  to  the 
cartesian  coordinates  {x,y,z)  are 

(39) 

(40) 

=  7K^f  - 
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(42) 


^^3) 

’'z  =  -  ^f>'{)  (44) 

=  jiy^^r,  -  y.z^)  (45) 

^y  = 

fz  =  -  ^>’{) 
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3  Numerical  Considerations 


Finite-Volume  Discretization 


With  the  aid  of  Equations  30-47,  all  derivatives  with  respect  to  x,  y,  and  z 
can  be  eliminated  from  the  governing  equations  in  Chapter  2,  leaving  only 
derivatives  with  respect  to  i;,  and  f.  For  discretization,  one  then  has  to 
define  the  grid  spacing  (A^,Ajj,Af)  in  the  curvilinear  (computational) 
coordinate  system.  The  latter  is  arbitrary,  and  standard  procedure  is  to 
impose  unit  spacing  for  convenience;  i.e.. 


A^  =  Ar,  =  Ar  =  1 


(48) 


Since  infinitesimal  volumes  in  the  two  coordinate  systems  are  related  by 


J  drj  =  dxdy  dz 


(49) 


it  then  follows  that  the  finite  volume  of  a  single  discrete  grid  cell  is  given  by 


J  = 


dx  dy  dz 


(50) 


In  computational  (^,Tj,f)  space,  all  the  grid  cells  are  unit  cubes;  but  in  car¬ 
tesian  {x,y,z)  space,  they  are  hexahedrons  of  arbitrary  shape,  as  shown  in 
Figure  1.  Thus,  a  triple  integral  appears  on  the  right-hand  side  of  Equa¬ 
tion  50,  with  the  limits  of  integration  being  the  six  faces  of  the  grid  cell.  All 
discrete  calculations  are  done  in  the  computational  (^.Jj.f)  space  for 
convenience. 

MAC3D  uses  grids  of  the  marker-and-cell  (MAC)  type  for  discretizing  the 
governing  equations  discussed  in  Chapter  2.  This  means  that  scalar  quantities 
such  as  pressure,  temperature,  density,  etc.,  are  defined  at  the  centers  of  the 
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grid  cells,  while  the  volumetric  fluxes  {U,V,W)  given  by  Equations  36-38  are 
defined  at  the  centroids  of  the  cell  faces.  Since  the  fluxes  are  proportional  to 
the  three  cartesian  components  of  velocity  the  latter  must  also  be 

defined  on  the  cell  faces.  In  practice,  one  needs  definitions  for  «,  v,  and  w  at 
the  cell  centers  as  well,  but  more  will  be  said  about  this  later. 

Discrete  flow  equations  are  obtained  for  each  cell  by  integrating  the 
governing  equations  over  the  cell  volume  and  converting  volume  integrals  to 
surface  integrals  via  the  Gauss  theorem.  The  details  of  this  procedure  are 
given  for  two  dimensions  in  Bernard  (1993),  which  discusses  the  depth- 
averaged  STREMR  model.  The  extension  to  three  dimensions  is  tedious  but 
straightforward. 


Computational  Coordinates 

Since  the  grid  has  unit  spacing  in  the  computational  space,  it  is 

convenient  to  replace  the  curvilinear  coordinates  with  the  integer 

computational  coordinates  (i,j,k).  Although  the  two  sets  of  coordinates  are 
identical,  the  integers  (i,j,k)  also  serve  as  indices  for  entire  grid  cells  and  cell 
faces,  in  addition  to  defining  the  locations  of  the  grid  nodes  in  computational 
space. 

In  the  MAC3D  indexing  scheme,  cell  is  the  particular  cell  whose 
corners  lie  inclusively  between  nodes  (M  j-1,^-1)  and  as  shown  in 

Figure  2.  For  example,  cefl  (1,1,1)  has  its  eight  corners  at  nodes  (0,0,0), 
(0,0,1),  (0,1,0),  (1,0,0),  (0,1,1),  (1,0,1),  (1,1,0),  and  (1,1,1). 

In  contrast  to  the  cells  themselves,  there  are  three  orientations  of  cell  faces. 
Faces  of  constant  i  are  called  i-faces:  faces  of  constant  j  are  called  j-faces;  and 
faces  of  constant  k  are  called  k-faces.  Face  indices  consist  of  one  node 
coordinate,  corresponding  to  the  constant  coordinate  on  that  particular  face, 
and  t^  c^  indices,  corresponding  to  the  remaining  two  indices  for  that 
particular  cell. 

Thus,  the  two  i-faces  for  cell  ii,j,k)  are  labelled  as  i-face  (i,j,k)  and  i-face 
{i-lj,k).  The  four  corners  of  i-face  {i,j,k)  lie  at  nodes  (i,j,k),  {i,j-l,k), 
{i,j,k-l),  and  (i,j-].,k-l).  The  four  corners  of  i-face  (i-lj,k)  lie  at  nodes 
{i-lj,k),  ii-\J-l,k),  (i-lj,k-l),  and  (i-lj-l,k-l)- 

For  example,  the  two  i-faces  of  cell  (1,1,1)  are  i-face  (1,1,1)  and  i-face 
(0,1,1).  I-face  (1,1,1)  has  its  four  corners  at  nodes  (1,1,1),  (1,1,0),  (1,0,1), 
and  (1,0,0).  On  the  opposite  side  of  the  cell,  i-face  (0,1,1)  has  its  four 
corners  at  nodes  (0,1,1),  (0,1,0),  (0,0,1),  and  (0,0,0). 
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Face-Centered  Fluxes 


Expressed  in  terms  of  the  curvilinear  coordinates,  the  incompressible 
constraint  for  conservation  of  mass  (Equation  3)  becomes 


dU  ^  dV  ^  dW 
dr} 


(51) 


where  U,  V,  and  W  are  the  volumetric  fluxes  given  by  Equations  36-38.  For 
a  MAC  grid  cell  in  computational  (i,j,k)  space,  U  is  defined  for  i-faces,  V  for 
j-faces,  and  ly  for  k-faces.  Using  the  same  indices  for  the  fluxes  and  the 
faces  of  cell  the  discrete  expressions  for  the  partial  derivatives  in 
are 


U{iJ,k)  -  m-l,j,k)  (52) 


V(iJ,k)  -  V(iJ-l,k)  (53) 


W(iJ,k)  -  W(iJ,k-l)  (54) 


In  simpler  terms,  the  volumetric  flux  through  a  cell  face  is  just  the  area  of  the 
face  multiplied  by  the  velocity  normal  to  the  face.  Thus,  a  mass-conserving 
flow  can  be  established  for  cell  (i,j,k)  by  correctly  defining  a  single  normal 
component,  instead  of  three  cartesian  components,  for  the  velocity  on  each  of 
the  six  faces. 


Equation  51 


drj 


dW 


Cartesian  Velocities 

If  all  three  velocity  components  («,  v,  and  w)  ^  known  for  each  cell  face, 
then  it  is  straightforward  to  calculate  the  appropriate  volumetric  flux  (U,  V,  or 
W)  through  each  face.  When  solving  Equations  36-38  for  u,  v,  and  w, 
however,  one  finds  that  each  of  these  velocity  components  requires  three  flux 
components  (U,  V,  and  W)  that  have  to  be  taken  from  three  distinct  cell  faces. 
Thus,  if  the  flux  alone  is  known  for  each  face,  there  is  not  enough 
information  to  calculate  all  three  velocity  components  at  a  single  location.  To 
compute  u,  v,  and  w  on  a  given  cell  face,  it  is  necessary  to  use  flux  infor¬ 
mation  from  other  cell  faces.  This  creates  a  dilemma  concerning  the 
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unambiguous  definition  of  velocity  components  in  the  equations  for  conser¬ 
vation  of  mass  (3  and  51)  and  momentum  (4). 

The  problem  is  as  follows.  To  maintain  conservation  of  mass  as  required 
by  Equation  51,  the  fluxes  (U,V,W)  their  time  increments  (AU,AV,AW) 
have  to  be  uniquely  defined  on  the  cell  faces.  But  the  flux  increments 
(AU,AV,AW)  depend  on  the  velocity  increments  (A«,Av,Aw)  via  Equa¬ 
tions  36-38,  and  these  cannot  be  uniquely  defined  at  a  given  location  unless 
the  velocity  components  themselves  are  uniquely  defined.  In  the  end,  one  has 
to  choose  between  two  alternatives:  averaging  flux  components  in  space  m 
averaging  flux  increments  in  time.  MAC3D  takes  the  latter  approach, 
implemented  in  a  predictor-corrector  scheme. 


Predictor-Corrector  Scheme  for  Updating  Fluxes 


In  the  predictor  phase,  flux  components  (U,  V,  W)  from  three  of  the  six 
faces  are  used  to  define  the  required  velocity  components  (u,v.w),  which  are 
used  in  Equation  12  to  calculate  provisional  velocity  increments 
(Am',Av',Aw').  These  are  inserted  into  Equations  36  -  38  to  obtain  non-mass- 
conserving  flux  increments  {AU',AV’,AW^,  which  are  then  added  to  the 
existing  flux  increments  to  obtain  provisional  flux  components  (U'.V'.W). 

The  latter  are  inserted  in  the  curvilinear  version  of  the  Poisson  equation  for 
pressure,  i.e.. 


div  grad  p  = 


Po 

JAt 


dU'  ^  ^  ^  dW' 
dr] 


(55) 


After  Equation  55  is  solved  iteratively  for  the  pressure  p,  using  a  precondi¬ 
tioned  conjugate  scheme  developed  by  Kapitza  and  Eppel  (1987),  the  resulting 
pressure  gradient  is  used  to  compute  the  correct,  mass-conserving  flux  incre¬ 
ments  (AU^KaV^KaW^^)  for  the  predictor  phase. 

The  corrector  phase  is  similar  to  the  predictor  phase,  except  that  the  flux 
components  {U,V,W)  from  the  other  three  cell  faces  are  used  to  define  the 
required  velocity  components  (u,v,w).  The  rest  of  the  procedure  for  com¬ 
puting  mass-conserving  flux  increments  (AU^‘^^,AV^^Ka\V^'^^)  for  the  corrector 
phase  is  the  same  as  that  for  the  predictor  phase.  At  the  end  of  the  corrector 
phase,  however,  the  flux  increments  obtained  in  both  phases  are  averaged  to 
give  the  net  increments  (AU,AV,AW)  for  one  time-step;  e.g.. 
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(56) 


Af/  =  +  A[/<‘'>) 


Note  that,  given  six  faces  to  choose  from,  there  are  eight  possible  com¬ 
binations  of  three,  whose  flux  components  can  be  used  to  define  velocity 
components  in  the  predictor  phase.  To  avoid  directional  bias,  MAC3D  uses 
each  of  the  eight  possible  combinations  once  in  every  eight  time-steps. 


Predictor-Corrector  Scheme  for  Transport 
Equations 

With  the  exception  of  the  turbulence  quantities  k  and  e  (to  be  discussed 
later),  a  predictor-corrector  scheme  is  used  for  updating  all  variables  whose 
governing  equations  have  the  general  form  of  Equation  14.  Expressed  in 
curvilinear  coordinates,  the  latter  becomes 


The  right-hand  side  of  Equation  57  contains  first-  and  second-order  derivatives 
in  space,  while  the  left-hand  side  contains  first-order  derivatives  in  space  and 
time.  In  both  the  predictor  and  corrector  phases,  the  time  derivative  is 
approximated  with 


^  ^  rp*  -  Ip  (58) 

dt  At 


where  \p*  is  the  new  value  and  \p  is  the  old  value.  In  finite-difference  jargon, 
this  is  called  forward-time  differencing. 

For  the  first-order  space  derivatives,  the  predictor  and  corrector  phases 
alternately  use  forward-space  and  backward-space  differencing.  For  example, 
if  the  predictor  phase  uses 


-  rP(iJ,k)] 


(59) 


then  the  corrector  phase  must  follow  with 
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(60) 


=  U(i-lJ,k)['P(iJ>k)  -  Hi-hj,k)] 
oz 


and  vice  versa.  For  all  three  coordinate  directions  there  are  eight  possible 
combinations  of  forward-  and  backward-space  differences  that  can  be  used  in 
the  predictor  phase.  MAC3D  uses  each  of  these  combinations  once  in  every 
eight  time-steps  to  avoid  accumulating  directional  bias. 

The  alternate  use  of  forward  and  backward  differencing  for  first-order 
space  derivatives  is  called  MacCormack’s  method  (MacCormack  1969; 
Bernard  1993).  In  contrast,  central  differencing  is  used  for  second-order  deri¬ 
vatives  in  space  in  both  the  predictor  and  corrector  phases;  e.g., 

-  2rP(iJ,k)  +  Hi-lJ,k)  (61) 


For  a  single  time-step,  the  net  increment  A\p  is  half  the  sum  of  the  predictor 
increment  A\p^^  and  the  corrector  increment  A\p(^^,  i.e.. 


A^P  =  +  A#) 

2 


(62) 


MAC3D  uses  the  foregoing  differencing  scheme  not  only  for  updating 
passive  constituents,  but  also  for  computing  provisional  velocity  increments 
(Au',Av',Aw')  via  Equation  12.  Expressed  in  its  curvilinear  form,  the  advec- 
tive  term  in  the  x-component  of  the  momentum  equation  becomes 


u  -grad  u  = 


I 

J 


U—  +  V—  +  W— 

dri 


(63) 


with  similar  expressions  for  u  ■  grad  v  and  u  ■  grad  w  in  the  y-  and  z- 
components,  respectively. 

In  Equation  12,  the  cartesian  velocity  components  {u,v,w)  are  treated  as 
though  they  were  cell-centered  quantities,  even  though  they  are  computed 
from  flux  components  {U,V,W)  taken  from  three  different  faces.  This  artifice 
allows  the  provisional  velocity  increments  (Am',Av',Aw')  to  be  computed  in 
the  same  manner  as  the  increments  of  truly  cell-centered  quantities.  Direc¬ 
tional  bias  is  small  (though  not  completely  insignificant)  after  each  full  time- 
step,  because  the  cartesian  velocity  components  {u,v,w)  are  computed  with 
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flux  components  {U,V,W)  taken  from  opposing  cell  faces  in  the  predictor  and 
corrector  phases. 

As  an  alternative  to  the  MacCormack  scheme  for  discretizing  u-grad  \p, 
MAC3D  also  offers  two-point  upwind  differencing  in  both  the  predictor  and 
corrector  phases.  For  example,  if  >  0,  then  upwind  differencing 

gives 


=  U{i-\j,k)U{ij,k)  -  i{j-\j,k)]  (64) 

OK 

Otherwise,  if  U(i,j,k)  <  0,  then  upwind  differencing  gives 

=  mj,k)  [  1  J,k)  -  ^p(ij,k)  ]  (65) 

oK 


Although  formally  less  accurate  than  the  MacCormack  scheme,  the  upwind 
scheme  guarantees  that  changes  in  \l/  propagate  only  in  the  direction  of  flow. 
It  is  intended  for  use  mainly  in  situations  where  the  MacCormack  scheme 
proves  numerically  unstable  or  otherwise  unsatisfactory. 


Single-Step  Upwind  Scheme  for  Turbulence 

The  governing  equations  for  turbulence  (22  and  23)  are  very  sensitive  to 
the  spatial  discretization  used  for  the  total  derivatives  Dk/Dt  and  De/Dt,  and 
also  to  that  used  for  the  velocity  derivatives  in  the  rate  coefficient  F  defined 
by  Equation  24.  After  considerable  experimentation  with  the  turbulence  equa¬ 
tions,  the  MacCormack  predictor-corrector  scheme  was  abandoned  in  favor  of 
a  single-step  upwind  scheme  for  advancing  k  and  e  in  time.  In  the  latter 
scheme,  both  u-grad k  and  u-grad  e  are  discretized  with  upwind  differ¬ 
encing,  and  their  increments  (AA:,Ae)  are  computed  in  a  single  (predictor)  step. 

In  contrast,  the  velocity  derivatives  appearing  in  F  cannot  be  approximated 
with  one-sided  differencing,  nor  can  the  velocities  themselves  be  calculated  in 
the  manner  used  in  the  predictor-corrector  scheme  for  the  momentum  equa¬ 
tion.  In  this  case,  the  flux  components  (U,V,W)  are  averaged  for  opposing 
cell  faces  to  obtain  truly  cell-centered  approximations  for  the  velocity 
components  {u,v,w).  The  cell-centered  velocities  are  then  used  in  central- 
difference  approximations  for  the  velocity  derivatives  in  F;  e.g.. 
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(66) 


The  foregoing  measures  are  necessary  to  keep  the  variation  of  k,  e,  and  >» 
smooth  and  well-behaved.  If  the  MacCormack  scheme  is  used  for  k  and  e,  or 
if  other  approximations  are  imposed  for  the  velocity  derivatives  in  F,  the 
computed  flow  may  become  unpredictably  erratic. 


Inflow/Outflow  Boundaries 

The  bottom,  sides,  and  free  surface  of  a  reservoir  are  true  boundaries. 

That  is,  they  are  impermeable,  and  no  fluid  passes  through  them.  In  contrast, 
inflow/outflow  boundaries  are  not  true  boundaries  at  all,  but  are  mathematical 
conveniences  that  facilitate  flow  computation  on  a  grid  of  limited  extent. 

MAC3D  recognizes  two  t3^es  of  inflow/outflow  boundaries  (FLUX  and 
OPEN),  both  of  which  allow  flow  fluid  to  enter  or  leave  the  grid.  The 
volumetric  flux  components  (U,V,W)  are  fixed  and  unchanging  (with  time)  on 
cell  faces  designated  as  FLUX  faces.  In  addition,  cell-centered  quantities  such 
as  temperature,  turbulence  energy,  etc.,  are  assigned  fixed  values  just  outside 
the  grid  along  FLUX  boundaries. 

OPEN  boundaries  are  frontiers  beyond  which  the  flow  is  unknown  and 
through  which  it  passes  as  if  no  boundary  were  present  at  all.  In  principle, 
the  only  property  of  OPEN  boundaries  is  transparency,  and  they  should 
transmit  outward-moving  disturbances  without  producing  any  reflections.  To 
obtain  boundary  values  for  velocity  components  normal  to  OPEN  faces,  as 
well  as  external  values  for  cell-centered  quantities  just  outside  OPEN 
boundaries,  MAC3D  uses  a  discrete  radiation  condition  proposed  by  Orlanski 
(1976).  Subject  to  this  boundary  condition,  normal  velocities  on  OPEN  faces 
are  free  to  change  in  response  to  the  flow  just  inside  the  grid,  and  cell- 
centered  quantities  just  outside  the  grid  are  free  to  change  in  a  similar  manner. 

Let  \p  represent  a  face-centered  (or  cell-centered)  quantity  whose  value  is 
needed  on  (or  just  outside)  an  OPEN  i-face.  The  discrete  radiation  condition 
imposes  the  following  equation  for  the  radiation  of  ^  through  the  i-face: 


\l/,  +  c\p^  =  0 


(67) 


Both  the  time  derivative  and  the  space  derivative  can  be  calculated  from 
the  existing  flow  in  the  previous  time-step  (or  predictor  phase).  The  one 
remaining  unknown  in  Equation  67  is  the  propagation  speed  c,  given  by 
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(68) 


The  derivatives  and  are  calculated  via  central-difference  approximations 
using  previous  values  for  ^  and  Axl/  at  locations  just  inside  the  grid. 

If  the  sign  of  c  is  such  that  the  direction  of  propagation  would  be  inward 
rather  than  outward,  the  computed  value  of  c  is  automatically  replaced  by 
zero.  Now  with  the  propagation  speed  completely  defined,  the  new  boundary 
value  '  is  given  by 

xj/^  =  xl/  -  c At x[/^  (69) 


in  which  the  derivative  xp^  is  computed  with  a  one-sided  (forward  or  back¬ 
ward)  difference  approximation  in  the  ^-direction.  For  example,  if  i-face 
is  an  OPEN  face,  and  an  external  value  is  needed  for  xp'  in  cell 
(i+lj,k),  then 


xp\i+lj,k)  =  xP(i^^lJ,k) 

-  cAt[xP(i+lJ,k)  -  ^(ij,k)] 


(70) 
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Composite  Grids 


MAC3D  uses  composite  grids  (also  known  as  multiblock  grids)  made  up  of 
one  or  more  subdomains.  These  can  have  any  shape  and  orientation  in 
cartesian  {x,y,z)  space,  and  each  subdomain  occupies  a  distinct  rectangular 
block  of  computational  space.  The  subdomains,  or  grid  blocks,  must  be 
contiguous  in  cartesian  space,  but  their  placement  in  computational  space  is 
arbitrary  as  long  as  they  are  separated  with  buffer  zones  at  least  two 
units  thick.  The  buffer  zones  are  needed  to  store  updated  values  of  flow 
variables  for  neighboring  cells  adjacent  to  the  subdomain  interfaces  (in 
cartesian  space). 

For  convenience,  each  grid  block  has  its  own  local  (i,j,k)  coordinates. 

Once  the  position  of  a  block  has  been  specified  in  the  global  (i,j,k)  coordinate 
system,  all  subsequent  input  for  that  block  refers  to  its  local  (i,j,k) 
coordinates.  Moreover,  if  a  calculation  is  to  be  made  for  a  limited  portion  of 
a  reservoir,  then  only  the  grid  blocks  that  encompass  the  region  of  interest 
need  to  be  loaded  into  the  calculation. 

Composite  grids  facilitate  calculations  in  which  fine  grid  spacing  is  needed 
near  a  particular  structure  or  boundary  feature,  but  not  elsewhere.  For 
example,  the  wake  created  by  a  cylindrical  pier  in  a  rectangular  channel  is 
best  computed  with  a  fine  cylindrical  grid  close  to  the  pier,  and  a  coarse 
rectangular  grid  far  downstream. 


Grid  Generation 


Grid  generation  can  be  accomplished  witn  any  of  several  3-D  numerical 
grid  generators  that  are  now  commercially  available.  The  main  things  to  keep 
in  mind  are,  first,  that  the  grid  nodes  have  to  match  perfectly  (with  no  over¬ 
lap)  on  the  common  boundaries  of  neighboring  subdomains;  second,  that 
variations  in  grid  spacing  should  be  gradual  rather  than  abrupt;  and,  third,  that 
extreme  skewness  of  intersecting  grid  lines  should  be  avoided.  As  a  rule  of 
thumb,  the  grid  spacing  should  change  by  no  more  than  30  percent  from  one 
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cell  to  the  next,  and  the  angle  between  intersecting  grid  lines  should  lie 
between  45  and  135  deg. 


Inflow/Outflow  Data 


MAC3D  accepts  user-specified  inflow/outflow  data  along  FLUX 
boundaries  and  OPEN  boundaries.  Input  velocity  distributions  may  be 
uniform,  linear,  or  parabolic.  Temperature  and  other  passive  constituents  may 
be  assigned  a  single  value  along  a  particular  boundary  segment,  or  otherwise 
extrapolated  from  the  initial  vertical  distribution  in  the  reservoir.  Uniform 
values  are  imposed  for  turbulence  energy  and  dissipation  rate  (or  eddy 
viscosity)  along  all  inflow/outflow  boundaries.  As  discussed  in  Chapter  3, 
initial  values  are  held  fixed  for  all  time  along  FLUX  boundaries,  but  they  are 
allowed  to  change  subject  to  the  discrete  radiation  condition  (e.g.. 

Equation  67)  along  OPEN  boundaries. 


Initial  Conditions 

MAC3D  uses  the  velocities  specified  along  inflow/outflow  boundaries  to 
compute  an  irrotational,  mass-conserving  flow  (potential  flow)  inside  the  grid, 
and  this  provides  an  initial  condition  for  velocity.  The  net  flow  rate  can  be 
either  computed  directly  from  the  inflow/outflow  velocities  or  specified 
independently.  In  the  latter  case,  MAC3D  adjusts  the  magnitudes  (but  not  the 
relative  distributions)  of  the  inflow/outflow  velocities  to  match  the  specified 
flow  rate.  The  initial  temperature  is  a  user-specified  vertical  distribution, 
which  can  be  linear,  cubic,  or  cubic  spline. 


Flow  Development 

Using  the  initial  and  boundary  conditions  previously  discussed,  MAC3D 
employs  the  numerical  algorithms  discussed  in  Chapter  3  to  solve  the  govern¬ 
ing  equations  outlined  in  Chapter  2.  The  discrete  solution  is  advanced  step  by 
step  through  time  to  produce  the  developing  flow.  If  a  steady  state  is  possible 
for  the  flow  under  consideration,  the  model  should  converge  to  that  state  and 
hover  about  it  (with  small  deviations  from  one  time-step  to  the  next). 

Because  of  the  alternating  directions  us-;d  in  the  predictor-corrector 
scheme,  small  deviations  from  the  steady  state  tend  to  repeat  themselves  over 
cycles  of  eight  time-steps.  Steady  state  is  usually  at  hand  when  the  maximum 
and  minimum  values  of  the  flow  variables  remain  unchanged  for  fifty  or  more 
of  these  cycles.  In  any  case,  when  the  real  flow  has  no  steady  state,  the 
computed  flow  likewise  has  no  steady  state. 
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The  size  of  the  time-step  At  has  to  be  specified  by  the  user,  but  MAC3D 
computes  and  prints  two  reference  values  that  can  be  used  for  guidance.  The 
first  of  these  is  the  maximum  allowable  time  increment  At^  based  on 
numerical  stability  considerations  (Bernard  1993).  The  second  is  the  period 
which  is  the  inverse  of  the  Brunt-Vaisala  frequency  N.  The  quantity 
is  the  fundamental  period  of  oscillation  for  a  stratified  fluid  in  the  presence  of 
gravity,  given  by 


PQ 


(71) 


As  a  rule  of  thumb,  the  time-step  should  be  slightly  less  than  and  at 
least  ten  times  smaller  than  t^. 

If  so  directed  by  the  input,  MAC3D  will  compute  and  store  the  initial 
flow,  and  then  stop.  At  this  point,  the  user  can  examine  the  printed  output, 
which  includes  quantities  like  At^^  and  t^,  and  set  At  to  a  suitable  value. 

The  computation  can  then  be  restarted  and  continued  for  a  specified  number  of 
time-steps,  after  which  At  can  be  reset.  This  procedure  can  be  repeated  as 
many  times  as  necessary. 


Flow  Visualization 

At  the  end  of  each  MAC3D  run,  all  the  major  flow  variables  are  stored  in 
output  files  that  can  be  used  with  appropriate  flow-visualization  software  to 
create  vectors,  streamlines,  color /contour  maps,  and  other  types  of  plots  from 
the  computed  results.  Output  from  MAC3D  is  generated  in  the  well-known 
PLOT3D  format,  which  is  accepted  by  several  different  visualization  packages 
that  are  now  commercially  available.  Although  the  flow  calculations  them¬ 
selves  can  be  executed  on  workstations  or  even  personal  computers  (PC’s), 
they  are  much  more  quickly  done  on  supercomputers.  In  the  latter  case,  the 
resulting  output  files  usually  have  to  be  transferred  back  to  a  workstation  or 
PC  for  flow  visualization. 
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Provisional  Grid  for  McCook  Reservoir 

Figures  3  and  4  show  overhead  views  of  the  surface  and  bottom  grids, 
respectively,  for  a  provisional  configuration  of  McCook  Reservoir.  The  north 
arrow  in  Figure  3  is  not  true  north;  it  is  defined  only  for  convenience  of 
orientation.  The  cutout  in  the  bottom  (Figure  4)  represents  a  sump.  The  grid 
as  a  whole  is  divided  into  ten  subdomains  (blocks),  which  together  occupy  a 
129  X  44  X  5  region  in  computational  space.  Figure  5  shows  an 
overhead  view  of  the  block  edges  (thin  white  lines)  with  the  sump  grid  inset, 
and  Figure  6  shows  a  perspective  view  with  the  sump  walls  and  bottom 
shaded.  Figure  7  offers  a  similar  view  of  the  reservoir  bottom  with  the  sump 
excluded. 

The  simulated  reservoir  has  a  maximum  width  of  3,260  ft  (east  to  west  in 
Figure  3)  and  a  maximum  length  of  4,655  ft  (north  to  south  in  Figure  3). 
When  filled  to  capacity,  it  is  113  ft  deep  along  the  ridge  of  the  central  berm 
(Figure  7),  135  ft  deep  along  the  plateau  immediately  to  the  west  of  the  berm, 
and  150  ft  deep  in  the  trough  on  the  extreme  west  side  of  the  grid.  East  of 
the  berm,  the  reservoir  is  uniformly  130  ft  deep  except  in  the  sump,  where  the 
total  depth  is  230  ft.  The  ramp  at  the  north  end  of  the  berm  (Figure  7) 
reduces  the  depth  gradually  to  a  minimum  value  of  20  ft. 

Each  grid  block  is  five  cells  deep,  including  the  block  that  represents  the 
sump.  Thus,  the  discrete  flow  field  is  five  cells  deep  except  for  the  region 
containing  the  sump  (Figure  8).  This  region  consists  of  two  blocks  stacked 
vertically,  and  it  is  ten  cells  deep  altogether.  In  the  MAC3D  flow  simula¬ 
tions,  as  discussed  in  the  following  section,  flow  was  injected  uniformly 
through  a  port  consisting  of  the  bottom  two  rows  of  cell  faces  on  the  right 
side  of  the  sump  (Figure  9),  and  withdrawn  uniformly  through  the  reservoir 
surface  (Figure  3).  The  inflow  port  and  the  reservoir  surface  were  designated 
as  FLUX  boundaries  in  the  MAC3D  input. 

The  computational  grid  is  fine  enough  to  capture  the  largest  eddies  in  the 
flow,  but  not  fine  enough  to  achieve  grid  independence  (i.e.,  negligible 
change  in  the  computed  flow  with  further  grid  refinement).  For  the  latter 
purpose,  ten  to  twenty  additional  grid  spaces  would  be  needed  in  the  vertical 
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for  each  subdomain.  In  some  regions  of  the  reservoir,  a  comparable  increase 
in  the  horizontal  resolution  might  also  be  necessary.  For  qualitative  model 
demonstrations,  however,  the  grid  used  here  is  adequate. 


Imposed  Flow  Conditions 

The  reservoir  was  assumed  to  be  completely  full.  The  velocity  through  the 
inflow  port  was  fixed  at  1  fps,  which  created  an  inflow  of  10,870  cfs.  This 
was  balanced  by  an  equal  outflow  with  a  uniform  vertical  velocity  of 
0.001  fps  at  the  reservoir  surface.  The  walls  and  bottom  of  the  reservoir 
were  given  a  Manning  coefficient  of  0.02,  and  the  turbulence  model  was 
activated  in  all  cases.  The  inflow  turbulence  energy  was  2  percent  of  the  total 
inflow  energy,  and  the  kinematic  eddy  viscosity  was  set  at  2  ft^/sec  along  the 
inflow  port.  The  Prandtl  number  was  set  to  unity,  making  the  thermal  dif- 
fusivity  equal  to  the  eddy  viscosity.  MAC3D  flow  simulations  were  then 
made  for  the  following  conditions: 

a.  Cold  inflow  (10  °C)  into  cold  reservoir  (10  °C) 

b.  Cold  inflow  (10  “C)  into  warm  reservoir  (20  °C) 

c.  Cold  inflow  (10  °C)  into  stratified  reservoir  (10  to  20  °C) 

(1)  With  background  heat  conduction 

(2)  Without  background  heat  conduction 

For  the  stratified  reservoir,  the  initial  temperature  profile  (Figure  10)  was 
20  °C  from  the  surface  to  a  depth  of  39  ft;  linear  with  depth  from  39  ft  to 
91  ft;  and  10  °C  below  91  ft. 

The  initial  flow  used  for  all  cases  was  potential  flow.  This  was  computed 
automatically  by  MAC3D  to  establish  an  irrotational,  mass-conserving  velocity 
field  throughout  the  entire  reservoir.  The  governing  equations  were  then 
marched  through  time  for  7,200  sec  (2  hours)  in  increments  of  1  sec  per  time- 
step.  By  this  time,  the  cold  reservoir  with  cold  inflow  (neutrally  buoyant 
flow)  had  essentially  reached  steady  state.  The  other  (buoyant)  cases  did  not 
reach  steady  state,  partly  because  cold  water  was  continuously  flowing  into  the 
reservoir  and  displacing  warmer  water.  The  stratified  case  was  run  with  and 
without  background  heat  conduction  (as  discussed  in  Chapter  2)  to  illustrate 
the  extent  to  which  the  adiabatic  boundary  conditions  and  thermal  diffusivity 
can  affect  the  computed  flow. 
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Presentation  of  Computed  Results 


To  facilitate  interpretation  and  qualitative  comparison  of  the  computed 
results,  the  same  information  is  presented  graphically  for  each  set  of  imposed 
flow  conditions.  First  comes  a  perspective  view  of  selected  velocity  profiles 
over  the  entire  reservoir  (e.g.,  Figure  11).  Next  are  shown  close-ups  of  the 
velocities  in  and  above  the  sump  (e.g.,  Figures  12  and  13),  followed  by  an 
overhead  view  of  the  surface  velocities  (e.g..  Figure  14).  Then  comes  a 
series  of  overhead  views  of  flow  patterns  in  five  successive  grid  planes  (e.g.. 
Figures  15-19),  beginning  at  the  surface  and  ending  at  80  percent  of  the  total 
depth  (to  the  top  of  the  sump).  These  flow  patterns  were  obtained  by 
normalizing  all  the  velocity  vectors  so  that  each  has  a  magnitude  of  unity. 
Thus,  the  displayed  vectors  indicate  only  the  local  direction  of  the  flow.  In  an 
overhead  view  for  a  single  grid  plane,  however,  the  projected  vector  lengths 
will  not  all  be  the  same  unless  the  vertical  velocity  component  is  uniform. 


Cold  Reservoir  with  Cold  Inflow 

Most  of  the  action  for  the  cold  reservoir  (Figures  11-19)  is  confined  to  the 
region  near  the  sump.  Strong  eddies  are  evident  in  both  the  elevation  and 
horizontal  planes  of  the  sump  itself.  These  are  complemented  by  a  strong 
eddy  just  to  the  south,  and  a  weaker  eddy  (near  the  bottom)  just  to  the  north. 
Velocities  in  the  southern  half  of  the  reservoir  are  very  small  in  comparison 
with  those  in  and  around  the  sump. 


Warm  Reservoir  with  Cold  Inflow 


Results  for  the  warm  reservoir  (Figures  20-28)  are  dramatically  different 
from  those  for  the  cold  (neutrally  buoyant)  reservoir.  In  this  case,  the 
downward  buoyant  force  increases  the  energy  required  for  the  cold  inflow  to 
penetrate  vertically  the  warmer  (lighter)  water  in  the  reservoir.  Thus,  the 
colder  the  inflow,  the  greater  the  tendency  for  water  leaving  the  sump  to  flow 
under  the  ambient  water  (as  in  Figure  21)  instead  of  up  and  through  it  (as  in 
Figure  12). 

The  underflow  away  from  the  sump  along  the  bottom  induces  a  return  flow 
along  the  surface,  as  shown  in  Figures  21  and  22.  At  the  elapsed  time 
(2  hours)  for  which  results  are  shown,  the  highest  velocities  are  confined 
roughly  to  the  northern  third  of  the  reservoir.  Questions  of  magnitude  aside, 
the  overall  flow  patterns  are  more  complex  than  those  for  the  neutrally 
buoyant  flow,  with  each  grid  plane  exhibiting  horizontal  recirculation  in  dif¬ 
ferent  locations  (Figures  24-28). 
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Stratified  Reservoir  with  Cold  Inflow 


Distinct  simulations  were  made  for  the  stratified  reservoir  with  and  without 
background  heat  conduction  (Figures  29-46).  These  should  be  regarded  as 
bounding  calculations  rather  than  as  realistic  simulations.  They  are  presented 
mainly  to  demonstrate  the  disparity  of  prediction  that  can  arise  because  of 
uncertainty  concerning  the  internal  rate  of  heat  transfer. 

Near  the  sump,  the  velocities  for  the  stratified  reservoir  with  background 
heat  conduction  (Figures  29-32)  loosely  resemble  those  for  the  unstratified 
warm  reservoir  (Figures  20-23).  In  this  case,  however,  there  is  a  strong 
vertical  recirculation  west  of  the  central  berm,  which  was  not  evident  for  the 
unstratified  reservoir.  The  flow  in  this  region  is  generally  east  to  west  along 
the  surface  and  west  to  east  along  the  bottom. 

The  flow  in  the  stratified  reservoir  without  background  heat  conduction 
(Figures  38-46)  is  totally  different  from  that  for  the  reservoir  with  background 
heat  conduction  (Figures  29-37).  Here  the  largest  velocities  are  confined  to 
the  region  near  the  sump  (Figures  38-41),  with  multiple  zones  of  horizontal 
recirculation  occurring  in  different  planes  at  various  distances  from  the  sump 
(Figures  42-46). 

The  difference  between  these  two  flows  indicates  the  relative  influence  of 
the  eddy  diffusivity.  If  the  initial  diffusivity  of  2  ft^/sec  were  held  constant,  it 
would  be  sufficient  to  bring  the  entire  reservoir  to  thermal  equilibrium  within 
2  hours  (in  the  absence  of  any  flow).  Thus,  the  flow  simulation  with  back¬ 
ground  heat  conduction  combined  the  effects  of  cold  inflow  and  rapid  internal 
heat  conduction.  Although  the  reservoir  was  initially  stratified,  it  gradually 
became  unstratified  far  from  the  sump,  with  a  temperature  of  about  15  °C. 

In  the  flow  simulation  without  background  heat  conduction,  the  local 
temperature  changed  only  in  response  to  the  inflow,  which  was  too  weak  to 
destratify  the  entire  reservoir  through  fluid  motion  alone.  In  this  case,  the 
buoyant  forces  were  strong  enough  to  restrict  the  main  action  to  the  region 
near  the  sump,  and  the  initial  temperature  distribution  hardly  changed  at  all. 

Figure  47  shows  a  gray-scale  map  of  the  initial  stratified  temperature  (in 
degrees  centigrade)  on  vertical  (elevation)  surfaces  inside  the  reservoir  and 
along  its  sides.  Figure  48  shows  the  temperature  map  for  the  flow  without 
background  heat  conduction  after  2  hours  elapsed  time,  and  Figure  49  shows 
the  corresponding  map  for  the  flow  with  background  heat  conduction.  For 
comparison.  Figure  50  shows  the  temperature  map  for  the  previous  case  of  the 
warm  reservoir  (20  °C)  with  cold  inflow  (10  “C)  after  2  hours  elapsed  time. 
Note  that  vertical  distances  have  been  magnified  by  a  factor  of  five  in 
Figures  47-50. 
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6  Conclusions  and 
Recommendations 


In  most  respects,  one  can  regard  MAC3D  as  a  three-dimensional  extension 
of  the  depth-averaged  STREMR  model.  MAC3D  uses  the  same  numerical 
algorithms  and  turbulence  model  as  STREMR,  but  it  also  includes  buoyancy 
and  offers  the  added  flexibility  of  composite  (multiblock)  grids.  In  any  case, 
MAC3D  draws  heavily  on  the  insight  and  experience  previously  gained  in 
developing  STREMR,  which  has  been  extensively  tested  and  documented 
(Bernard  1992,  1993;  Bernard  and  Schneider  1992;  Maynord,  in  preparation). 

This  report  offers  no  comparison  of  model  predictions  with  physical  test 
data,  but  MAC3D  is  not  completely  unverified.  At  each  stage  of  develop¬ 
ment,  the  model  has  been  run  and  checked  against  STREMR  predictions  for 
2-D  flow.  Thus,  in  certain  3-D  settings  for  which  the  flow  should  be  2-D, 
MAC3D  actually  produces  a  2-D  flow;  and  the  predictions  agree  with  those  of 
STREMR.  In  future  work,  MAC3D  predictions  should  be  tested  against 
pertinent  experimental  results  for  fully  3-D  flow. 

The  existing  model  now  seems  ready  for  implementation  as  a  supplemen¬ 
tary  aid  in  design  studies  of  McCook  Reservoir,  although  undiscovered 
deficiencies  may  still  have  to  be  dealt  with.  Given  a  specified  inflow  or  out¬ 
flow,  MAC3D  can  predict  reasonable  bounds  for  the  resulting  circulation 
inside  the  reservoir.  The  major  uncertainties  are  the  rates  of  surface  heat 
exchange  and  internal  heat  conduction,  and  these  merit  considerable  further 
study.  In  particular,  modules  need  to  be  developed  to  account  for  reflection 
and  absorption  of  solar  radiation,  shear  stress  produced  by  surface  winds,  and 
surface  heat  flux  produced  by  changes  in  ambient  atmospheric  temperature. 

MAC3D  is  capable  of  predicting  the  hydrodynamic  response  to  inflow, 
outflow,  and  hydraulic  mixers,  but  ic  also  needs  to  predict  the  response  to 
bubble  diffusers.  The  latter  represents  a  special  problem,  because  rising  bub¬ 
bles  generally  move  with  a  different  velocity  from  that  of  the  flow.  Thus, 
future  work  should  include  the  development  of  adequate  models  for  localized 
bubble  plumes  and  their  influence  on  the  flow.  These  might  include  special 
buoyant  forces  added  along  the  vertical  planes  or  axes  of  diffusers,  along  with 
empirical  rules  for  the  spreading  of  the  plumes  themselves. 
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At  present  MAC3D  deals  with  hydrodynamics,  but  not  with  water  quality. 
To  make  the  model  more  generally  useful  as  an  aid  for  reservoir  design, 
operation,  and  modification,  future  work  should  include  the  addition  of 
modules  to  account  for  oxygen  transport  and  degradation,  as  well  as  pertinent 
chemical  and  biological  processes. 
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Figure  5.  Overhead  view  of  subdomain  edges  with  sump  grid  inset 
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Figure  8.  Cutaway  of  bottom  and  side  grids  in  and  around  sump 
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Figure  12.  Velocities  in  elevation  plane  of  sump  for  cold  reservoir  with  cold  inflow 
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Figure  15.  Surface  flow  patterns  for  cold  reservoir  with  cold  inflow 


Figure  17.  Flow  patterns  at  40  percent  of  total  depth  for  cold  resen/oir  with  cold  inflow 


Figure  18.  Flow  patterns  at  60  percent  of  total  depth  for  cold  reservoir  with  cold  inflow 
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Figure  19.  Flow  patterns  at  80  percent  of  total  depth  for  cold  reservoir  with  cold  inflow 


Figure  20.  Velocity  profiles  for  warm  reservoir  with  cold  inflow 


Figure  22.  Velocities  for  horizontal  planes  in  and  above  sump  for  warm  reservoir  with  cold 
inflow 


Figure  23.  Surface  velocities  for  warm  reservoir  with  cold  inflow 


Figure  25.  Flow  patterns  at  20  percent  of  total  depth  for  warm  reservoir  with  cold  inflow 


Figure  26.  Flow  patterns  at  40  percent  of  total  depth  for  warm  reservoir  with  cold  inflow 
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Figure  27.  Flow  patterns  at  60  percent  of  total  depth  for  warm  reservoir  with  cold  inflow 


Figure  28.  Flow  patterns  at  80  percent  of  total  depth  for  warm  reservoir  with  cold  inflow 


Figure  29.  Velocity  profiles  for  stratified  reservoir  with  cold  inflow 


Figure  30.  Velocities  in  elevation  plane  of  sump  for  stratified  reservoir  with  cold  inflow 


Figure  31 .  Velocities  for  horizontal  planes  in  and  above  sump  for  stratified  reservoir  with 
cold  inflow 
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Figure  33.  Surface  flow  patterns  for  stratified  reservoir  with  cold  inflow 


Figure  34.  Flow  patterns  at  20  percent  of  total  depth  for  stratified  reservoir  with  cold 
inflow 
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Figure  35.  Flow  patterns  at  40  percent  of  total  depth  for  stratified  reservoir  with  cold 
inflow 


"rCx  s  \  i  / 

'^^1''"  ">’i> f' '/ jljf/imK 

'^mm 

/iiifi  !  '" ' " '" ! ' "  •'« : 


Mmm 


V  ^  /  y  ^  /y 

" C  ^ ^  v^' ' \ J [ ^  ^ ^ '  /  / / < 


,^'^VNNVS- - 1 

>///H  Wnvv..^. 
z:^////  V'^^Nvy. 


Figure  37.  Flow  patterns  at  80  percent  of  total  depth  for  stratified  reservoir  with  cold 
inflow 


Figure  39.  Velocities  in  elevation  plane  of  sump  for  stratified  reservoir  with  cold  inflow  and 
no  background  heat  conduction 


Figure  40.  Velocities  for  horizontal  planes  in  and  above  sump  for  stratified  reservoir  with 
cold  inflow  and  no  background  heat  conduction 


isn'xsmsseasiiv-: 


li^liiSIliit® 


^?.-i;>>bri-ifV-:; 


Figure  41 .  Surface  velocities  for  stratified  reservoir  with  cold  inflow  and  no  background 
heat  conduction 


Figure  42.  Surface  flow  patterns  for  stratified  reservoir  with  cold  inflow  and  no  background 
heat  conduction 


Figure  43.  Flow  patterns  at  20  percent  of  total  depth  for  stratified  reservoir  with  cold 
inflow  and  no  background  heat  conduction 


Figure  44.  Flow  patterns  at  40  percent  of  total  depth  for  stratified  reservoir  with  cold 
inflow  and  no  background  heat  conduction 


Figure  45.  Flow  patterns  at  60  percent  of  total  depth  for  stratified  reservoir  with  cold 
inflow  and  no  background  heat  conduction 


Figure  47.  Gray-scale  map  of  initial  temperature  for  stratified  reservoir 


Figure  49.  Gray-scale  map  of  final  temperature  for  stratified  reservoir  with  background  heat 
conduction 


Appendix  A 
Notation 


D/Dt  Substantive  (total)  time-derivative  operator 

div  Divergence  operator 

g  Acceleration  due  to  gravity 

grad  Gradient  operator 

i,j,k  Integer  computational  coordinates 

J  Jacobian  of  coordinate  transformation 
k  Turbulence  energy 

n  Manning’s  coefficient 

n  Unit  vector  normal  to  boundary 

p  Pressure 

Pj  Ratio  of  thermal  diffiisivity  to  kinematic  viscosity,  also 
known  as  the  Prandtl  number 

S  Source/sink  term 

t  Time 

T  Temperature 

Tg  Static  (background)  temperature 

u  Vector  velocity 

u,v,w  Cartesian  components  of  velocity 
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U,V,W  Volumetric  fluxes  in  the  rj-,  and  f-directions,  respectively 
x,y,z  Cartesian  coordinates 
e  Turbulence  dissipation  rate 
6  Relative  deviation  from  reference  value 
IJL  Sum  of  the  molecular  viscosity  and  the  eddy  viscosity  iij- 
V  Kinematic  viscosity 
Diffiisivity  for  ^ 

Curvilinear  coordinates 
p  Density 

Pq  Reference  (maximum)  density 

T  Shear  stress 

rp  Scalar  function  or  vector  component 
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continuity  equation.  To  enhance  geometric  flexibility,  the  model  uses  composite  structured  grids  assembled  from 
one  or  more  curvilinear  subdomains.  Flow  simulations  are  presented  without  confirmation  for  a  provisional 
design  of  McCook  Reservoir.  Empirical  validation  of  the  model  will  be  undertaken  in  future  work. 
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